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1  Introduction 

Rotor  systems  are  often  constructed  using  either  rolling  element  or  journal  bearings  to  enable 
relative  angular  motion  between  the  rotor  shaft  and  the  system  housing.  Rotor  dynamics  and 
structural  finite  element  models  of  such  systems  generally  use  linear  elastic  elements  (linearized 
about  a  normal  operating  point)  to  represent  the  coupling  characteristics  of  these  components 
instead  of  employing  complicated  bearing  models  that  must  be  solved  simultaneously  with  the 
system  model.  The  process  of  determining  the  coupling  characteristics  of  both  types  of  bearings 
can  be  complicated,  but  for  different  reasons.  The  difficulty  for  rolling  element  bearings  is 
mainly  a  result  of  having  a  discrete  number  of  rolling  elements  (see  Figure  1.1a),  which  causes 
the  stiffness  to  be  highly  nonlinear  since  the  number  of  rolling  elements  in  contact  with  the 
raceway  at  any  given  time  depends  on  the  amount  of  shaft  deflection.  A  detailed  discussion  of 
rolling  element  bearing  modeling  is  out  of  the  scope  of  this  report.  Instead,  the  reader  is  referred 
to  the  discussion  of  this  topic  provided  in  Reference  1. 

Unlike  rolling  element  bearings,  journal  bearings  are  difficult  to  model  because  their  stiffness 
and  damping  depends  on  a  hydrodynamic  lubrication  layer  with  characteristics  that  are  largely 
controlled  by  the  relative  velocity  of  the  bearing  and  journal  surfaces,  lubrication  viscosity,  and 
bearing  geometry.  This  layer  exists  in  the  fluid  region  between  the  shaft  and  the  journal  as 
shown  in  Figure  1.1b.  The  ability  of  the  fluid  film  to  carry  a  load  is  a  result  of  a  converging  fluid 
film  gap  caused  by  eccentricity  of  the  journal  relative  to  the  bearing  and  the  presence  of  a 
viscous  fluid.  The  result  of  these  conditions  is  a  viscous  shear  flow  that  yields  a  net 
hydrodynamic  force  in  the  gap  that  can  support  a  load  with  no  direct  joumal-to-bearing  contact 
(and  hence  no  bearing  wear). 

There  are  a  variety  of  hydrodynamic  journal  bearing  designs  used  for  rotating  equipment  [2], 
most  of  which  are  a  variation  of  the  common  plain  bearing  design.  Some  of  the  usual  variations 
include  the  multi-lobed  and  tilting  pad  designs  shown  in  Figure  1.2.  Often,  the  intent  of  the 
design  modifications  is  to  enhance  the  overall  stability  of  the  rotor  system.  From  a  modeling 
perspective,  the  modifications  result  in  a  more  complex  system.  For  low  order  (1 -dimensional) 
shaft/housing  models,  methods  are  available  to  derive  net  stiffness  and  damping  characteristics  at 
the  shaft  centerline.  Higher  order  models  that  attempt  to  simulate  the  dynamic  behavior  around 


(a) 


(b) 


Figure  1.1.  (a)  Rolling  element  bearing  schematic,  and  (b)  plain  journal  bearing  schematic. 
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the  housing  and  shaft  circumference  require  distributed  stiffness  and  damping  representations  of 
the  bearings.  For  example,  a  rim-drive  motor  finite  element  model  may  require  the  application 
of  circumferentially  varying  forces  or  the  prediction  of  response  harmonics  in  the  circumferential 
direction  and  thus  a  distribution  of  the  bearing  coefficients  around  the  journal  circumference  is 
required.  It  is  the  intent  of  this  report  to  summarize  the  methods  currently  used  to  approximate 
the  equivalent  dynamic  characteristics  for  the  plain  journal  bearing  and  the  tilting  pad  bearing 
and  introduce  a  method  to  circumferentially  distribute  these  coefficients  for  use  in  a  structural 
finite  element  model.  As  previously  stated,  the  plain  journal  bearing  and  the  tilting  pad  bearing 
are  being  focused  on  here  because  they  are  two  of  the  most  commonly  employed  bearings  for 
rotating  equipment. 


Lobes 


Journal 


raas 


Housing 


Journal 


Pad  Pivot 


(a)  (b) 

Figure  1.2.  (a)  Four-lobed  bearing,  and  (b)  tilting  pad  bearing. 

2  Background 

The  process  of  computing  bearing  dynamic  coefficients  begins  with  solving  for  the  pressure  field 
in  the  fluid  film  between  the  journal  and  bearing  surfaces.  The  equation  used  for  this  purpose  is 
referred  to  as  the  Reynolds  equation.  Important  aspects  of  this  partial  differential  equation  are 
highlighted  next,  along  with  boundary  conditions  and  methods  to  solve  the  system. 

2.1  The  Reynolds  Equation 

The  film  pressure  in  a  journal  bearing  is  generally  computed  using  the  Reynolds  equation  shown 
below1  [4,5,6, 7], 


1  An  informative  historical  review  of  hydrodynamic  lubrication  theory  and  the  development  of  the  Reynolds 
equation  can  be  found  in  Reference  3. 
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where  h  is  the  film  thickness,  p  is  the  lubricant  viscosity,  p  is  the  fluid  film  pressure,  R  is  the 
bearing  radius,  t  is  time,  Q  is  the  angular  speed  of  the  journal  relative  to  the  bearing,  z  is  the 
coordinate  along  the  axis  of  rotation,  and  6  represents  rotation  about  z.  The  derivation  of 
Equation  2.1  begins  with  the  Navier-Stokes  and  continuity  equations  and  employs  several 
assumptions,  including  the  following:2 

•  effects  of  bearing  curvature  are  negligible  (i.e.,  centrifugal  forces  are  neglected), 

•  pressure  variation  across  the  fluid  film  is  negligible, 

•  the  flow  is  laminar, 

•  fluid  inertia  effects  are  negligible  when  compared  to  the  viscous  effects,3 * 

•  the  flow  is  incompressible,  and 

•  the  fluid  viscosity  is  constant. 


For  turbulent  flow,  turbulent  correction  factors,  Gg  and  Gz,  obtained  from  Hirs  bulk-flow  theory 
[9,10]  are  often  added  to  Reynolds  equation  to  approximately  account  for  the  turbulence  effects. 
Using  these  factors,  the  Reynolds  equation  becomes: 
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(2.3) 

(2.4) 


For  flows  with  smooth  surfaces  and  Reynolds  numbers  (Re  =  -^^)  less  than  105,  common 

P 

values  for  n  and  m  are:  n  =  0.066  and  m  =  -0.25  [11].  For  laminar  flow,  G<?and  Gz  are  both  equal 
1 

tol2- 


2  Inherent  to  the  Navier-Stokes  equations  is  the  assumption  of  a  Newtonian  fluid.  For  non-Newtonian  fluids,  a 
lubrication  rheological  model  may  be  necessary  (see  Reference  8  for  an  introduction  to  lubrication  rheology). 

3  The  consequence  of  this  assumption  is  that  all  terms  in  the  Navier-Stokes  equations  that  are  multiplied  by  f'Re  are 

neglected  when  ^Re  «  1  (£  =  HIL ,  Re  =  />V*H //4  H  is  a  characteristic  film  thickness,  L  is  the  bearing  width,  V  is  a 

characteristic  velocity,  and  /i  is  a  characteristic  fluid  viscosity). 
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Note  that  the  derivation  of  the  Reynolds  equation  assumes  negligible  fluid  inertia  effects.  This 

assumption  is  supported  when  the  modified  Reynolds  number,  Re*  =  — ,  is  much  less  than 

H  R 

unity.  For  cases  when  this  criterion  is  not  met,  methods  presented  in  References  10  and  16 
should  be  used  to  account  for  the  inertia  effects.4  Also  note  that  the  Reynolds  equation  assumes 
constant  fluid  viscosity  within  the  fluid  film.  Methods  to  incorporate  a  variable  viscosity  due  to 
thermal  effects  are  presented  in  References  17,  18,  and  19. 


2.2  Boundary  Conditions 

The  choice  of  appropriate  boundary  conditions  for  the  Reynolds  equation  has  been  the  topic  of 
much  discussion  over  the  years,  focusing  mainly  on  how  the  film  rupture  boundary  condition  is 
handled.  In  fact,  boundary  conditions  that  are  compatible  to  both  the  boundary  pressure  and 
continuity  are  yet  to  be  defined  [16].  Film  rupture  may  occur  at  the  ‘downstream’  boundary  if  a 
diverging  film  exists  at  this  location.  For  a  plain  journal  bearing,  a  diverging  film  will  always 
occur  at  the  downstream  boundary  (02=  n+<p  \r\  Figure  2.1).  However,  short-length  pads  often 
have  a  converging  film  over  the  entire  surface  (see  Figure  2.2  for  a  schematic  of  a  bearing  with 
short  length  pads). 


There  are  two  types  of  boundary  conditions  that  are  commonly  used:  the  Reynolds  conditions 


r=  o 


e:  Eccentricity 
O.  Journal  Diameter 
h:  Film  Thickness 
y  W.  Applied  Load 
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£2.  Journal  Speed 
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Figure  2.1.  Plain  journal  bearing  variables. 


4  Reference  10  reports  that  as  CIR  decreases,  the  results  from  the  equations  that  include  inertia  effects  decrease 
monotonically  to  the  results  of  the  classical  theory.  Thus,  it  may  be  acceptable  to  neglect  inertia  effects  when  CIR  is 
small. 
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Applied  Load,  W 


Figure  2.2.  Tilting  pad  journal  bearing  analysis  variables. 


dp_ 

dd 


=  0)  and  Gumbel’s  conditions. 


The  Reynolds  conditions  are  applied  by  setting 


every  term  of  negative  pressure  to  zero  during  the  solution  process,  and  Gumbel’s  conditions 
suggest  that  all  pressures  in  the  solution  that  are  less  than  atmospheric  pressure  be  ignored5  [16]. 
While  neither  method  maintains  flow  continuity,  the  theoretical  results  obtained  from  using 
either  of  these  conditions  have  been  shown  to  closely  match  experimental  data  [16]. 


2.3  Reynolds  Equation  Solution 

Closed  solutions  to  the  Reynolds  equation  (Equation  2.1)  have  been  developed  for  two  simplified 
conditions  [4]:  an  infinitely  long  bearing  ( Sommerfeld  or  long  bearing  solution),  and  an  infinitely 
short  bearing  ( Ocvirk  or  short  bearing  solution),  both  of  which  are  summarized  in  the  Appendix. 
The  Sommerfeld  solution  is  applicable  to  bearings  with  length  (L)  to  diameter  (D)  ratios  large 
enough  that  leakage  from  the  bearing  ends  is  small  relative  to  the  flow  in  the  circumferential 
direction  (i.e.,  UD  »  1),  whereas  the  Ocvirk  solution  is  applicable  to  bearings  with  LID  <  ~1 
such  that  end  leakage  is  significant.  Because  most  bearings  in  use  have  UD  ratios  in  the  range  of 
14  to  1,  the  Ocvirk  solution  is  often  applicable  [4]. 

In  order  to  avoid  the  limiting  assumptions  of  these  simplified  solutions,  the  Reynolds  equation  is 
often  solved  numerically  using  finite  difference  [12,13],  finite  element  [11],  or  other  methods 
wherein  the  solution  is  obtained  from  a  search  for  the  minimum  value  of  the  functional  [11]: 


5  This  is  sometimes  called  the  half-Sommerfeld  approach. 
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Moreover,  the  computation  is  sometimes  simplified  by  using  “Shelly’s  Method”,  which  assumes 
a  functional  dependence  of  the  pressure  in  the  axial  direction  [1 1,14]: 


p(0.z)=P,  1-y 


(2.6) 


Here,  n  is  a  chosen  constant  which  depends  on  the  length-to-diameter  ratio  of  the  bearing6 7,  and 


3  Fixed  Pad  Dynamic  Coefficients 

The  calculation  of  dynamic  coefficients  for  plain  and  fixed  pad  journal  bearings  is  considered  in 
this  section.  The  calculation  for  the  tilting  pad  bearing  builds  upon  the  fixed  pad  method,  as 
presented  in  Section  4. 

3.1  Journal  Load  Capacity 

Once  the  pressure  distribution  is  determined  using  one  of  the  methods  described  in  Section  2,  the 
external  load  that  can  be  supported  by  the  journal  is  determined  by  integrating  the  pressure 
distribution  over  the  film  surface  [5,13,16]: 


(3.1) 


where  Fx  is  the  load  in  the  jc-direction,  and  Fy  is  the  load  in  the  y-direction  (see  the  coordinate 
system  in  Figure  2.1  or  Figure  2.2).  The  integration  limits  for  Equation  3.1  correspond  to  the 
beginning  of  (#/)  and  the  end  of  ( 02 )  the  active  film  region  since  the  pressures  are  zero  elsewhere 
when  the  boundary  conditions  described  in  Section  2.2  are  applied/  Integrating  over  the  entire 
film  yields  the  net  load  that  can  be  supported  by  the  journal  at  the  specified  operating  conditions 
for  which  the  pressure  field  was  computed.  The  force  due  to  a  portion  of  the  overall  film  is 
computed  in  a  similar  manner,  but  with  integration  limits  specific  to  the  region  of  interest. 

3.2  Perturbation  Method  for  Computing  Dynamic  Coefficients 

The  reaction  forces  computed  in  Equation  3.1  can  be  related  to  dynamic  stiffness  and  damping 
coefficients  by  expanding  these  forces  about  the  static  offset  position  (i.e.,  the  normal  operating 
point)  using  a  Taylor’s  series  expansion  [5,13]: 


6  For  medium-length  journal  bearings  where  V*  <  UD  <  1,  n  =  2  provides  reasonable  accuracy  [  1 1,14). 

7  di  and  0:  are  limited  to  0<  9<  (<?+n)  for  a  plain  journal  bearing  based  on  the  converging  film  region. 


Fx  =  Fx0  +  kM  +  k^Ay  +  b^Ax  +  b^Ay 

Fy  =  Fy0  +  M*  +  kyyty  +  KM  +  byyAy  ’ 


(3.2) 


where  F ^  is  the  jr-component  of  the  static  load,  Fyo  is  the  ^-component  of  the  static  load.  Ax  and 
Ay  are  small  amplitude  motions  in  the  x-  and  y-directions,  respectively,  and  Ax  and  Ay  are  the 
vibration  velocities  in  the  x-  and  y-directions,  respectively.  The  coefficients  of  the  Ax,  Ay,  Ax , 
and  Ay  terms  in  Equation  3.2  constitute  the  equivalent  bearing  dynamic  coefficients  used  for 
dynamic  analyses  [5].  These  coefficients  are  determined  by  solving  the  Reynolds  equation  using 
a  perturbation  on  the  film  thickness,  Ah,  such  that  [5,7]: 

h  =  h0  +  Ah,  (3.3) 


where 


ho  =  nominal  film  thickness  =  c  +  x0cos#+  y0-sin#. 

Ah  =  perturbed  film  thickness  =  zlx-cos(#)  +  Ay-sin(#),  and 


ch 

dt 


Ajccos#  +  Aysin#. 


(3.4) 

(3.5) 

(3.6) 


The  perturbed  film  thickness  creates  a  perturbed  film  pressure  that  is  defined  using  a  Taylor’s 
series  expansion  [5,8,16,23]: 

P  =  Po  +  PxAx  +  pyAy  +  p'xAx  +  p'yAy  =  po  +  Ap,  (3.7) 


where  px  and  py  are  the  pressure  components  related  to  the  film  stiffness  and  p'x  and  p'y  are  the 
pressure  components  related  to  the  film  damping. 


Substituting  Equations  3.3  to  3.7  into  Equation  2.1  and  retaining  first  order  terms  only  yields  five 
equations  (coefficients  of  Ax,  Ay,  Ax ,  and  Ay  are  equated  in  this  step): 


Z{p,}  =  -Q— 
1  0)  2  dd 


Z{pJ=-  3 

z{p,h- 3 


cos#1^3/i0 
h0  2  dO 

™°Lq?!Sl- 3. 

h0  2  d0 


K  *  fyo  1  d 


^  cos  0^ 


\2/u  R  dO  Rd0{ 

K  1  dp0  1  d  ( sin# 


y 


12 fi  Rd6  RdG 


\  ffo  y 


— —  £2sin# 
2 

-  — Qcos# 
2 


Z{/?'}=cos# 
z{py}=  sin# 


(3.8) 


where  the  left  hand  operator  Z{ }  is  defined  as: 
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(3.9) 


Z{  }  =  — — 
W  Rdd 


hi  i  a 


\ 


12//  R  d6 

By  combining  Equations  3.1,  3.2,  and  3.7  we  have: 

** 

K 

kyy 

bXX 

K 
K 

hyy 


dz 


12//  dz 


[cos#j 
[sin#  j 

ii 

i 

N: 

[cos#l 
[sin#  j 

i 

ii 

cos#] 

sin#} 

ii 

i 

[cos#] 

[sin#} 

(3.10) 


Integration  over  the  entire  film  surface  in  Equation  3.10  yields  the  dynamic  coefficients  that  are 
often  used  in  rotor  dynamic  analyses.  For  high  order  structural  finite  element  analyses,  these 
coefficients  can  be  distributed  along  the  journal  circumference  by  specifying  appropriate 
integration  limits  corresponding  to  the  finite  element  node  locations.  It  must  be  noted,  however, 
that  the  distributed  bearing  coefficients  do  not  account  for  changes  in  the  film  that  result  from 
flexure  of  either  the  journal  or  bearing,  and  thus  this  method  is  only  an  approximation  for  cases 
when  the  journal  and/or  bearing  are  flexible.  See  Reference  23  for  a  discussion  of  the  effect  of 
such  flexibility  on  the  dynamic  coefficients. 


4  Tilting  Pad  Dynamic  Coefficients 

Computation  of  dynamic  coefficients  for  the  tilting  pad  bearing  is  more  involved  than  what  is 
presented  in  the  previous  section  for  the  ‘fixed  pad’  bearings.  The  method  presented  below  for 
tilting  pad  bearings  builds  upon  the  methods  described  above. 

4.1  Tilting  Pad  Characteristics 

A  tilting  pad  bearing  is  composed  of  N  pads,  each  with  the  capability  of  pivoting  around  a 
fictitious  point  located  on  a  circle  called  the  pivot  circle.  The  main  parameters  include  (see 
Figure  2.2):  the  number  of  pads,  N;  the  journal  radius,  /?;  the  pivot  point  circle  clearance,  C';  the 
pad  angle,  y,  equivalent  mass  of  the  pad,  M  =  I/R2\  pivot  location,  pivot  pad  location  ratio, 
cdy\  pad  width,  L;  pad  radius  of  curvature,  Rp\  and  the  machined  clearance,  C  =  Rp-  R. 

The  preload8,  also  referred  to  as  the  preset,  is  an  important  analysis  parameter  for  the  tilting  pad 
bearing.  The  preload  is  a  measure  of  the  converging  film  that  is  created  by  the  bearing  geometry 
and  not  by  an  external  load  as  the  name  might  infer.  Extreme  bearing  preloads  are  shown  in 
Figure  4.1  to  aide  in  the  visualization  of  this  parameter. 


8  Note  that  the  preload,  m=  1  -  C'/C ,  is  different  from  the  preload  ratio,  C'IC. 


8 


Constant  Film 
Thickness 


Figure  4.1.  Concept  of  bearing  preload:  (a)  preload  =  1  (preload  ratio  =  0),  (b)  preload  =  0. 

An  important  characteristic  of  the  tilting  pad  bearing  that  separates  it  from  the  other  types  of 
bearings  is  the  existence  of  a  rotational  degree  of  freedom  for  each  pad.  These  additional 
degrees  of  freedom  result  in  a  total  of  4+5N  stiffness  and  4+5 N  damping  coefficients  [11]. 
Because  rotor  dynamic  models  generally  are  not  concerned  with  modeling  the  motion  of  the 
individual  pads,  methods  have  been  developed  to  obtain  a  reduced  set  of  coefficients  (4  stiffness 
and  4  damping)  by  assuming  a  vibration  frequency  for  the  pad  rotation.  A  discussion  of  reduced 
coefficients  is  provided  below  in  Section  4.3.  Prior  to  performing  any  type  of  dynamic 
coefficient  reduction,  however,  the  static  characteristics  must  be  defined. 

4.2  Pad  Assembly  Method 

The  procedure  for  computing  tilting  pad  dynamic  coefficients  begins  with  determining  static 
characteristics  of  each  tilting  pad  operating  as  a  fixed  pad.  The  characteristics  of  interest  are  the 
force  and  journal  eccentricity.  Once  the  static  characteristics  are  known  for  a  single  fixed  pad, 
the  static  characteristics  of  the  tilting  pad  bearing  are  determined  using  superposition  of  the  fixed 
pad  results.  This  method,  originally  introduced  by  Lund  [13]  and  used  by  other  researchers  in 
subsequent  studies  [20,21,22],  is  know  as  the  Pad  Assembly  Technique.  The  local  eccentricity, 
£„  and  attitude  angle,  #,  of  each  pad  is  determined  using  Equation  4.1. 

£ ,  cos <p,=m- e0  cos (y/t  -</>0),  (4.1) 

where  Ei  is  the  eccentricity  ratio  relative  to  pad  i,  $  is  the  attitude  angle  relative  to  pad  i,  is  the 
eccentricity  ratio  of  the  assembled  bearing,  and  fo  is  the  attitude  angle  of  the  assembled  bearing 
relative  to  the  global  coordinate  system  (see  Figure  4.2). 

Equation  4.1  is  derived  from  geometrical  considerations  using  a  local  pad  coordinate  system  as 
shown  in  Figure  4.2  [13].  This  assembly  method  generally  involves  a  trial  and  error  approach  to 
find  the  assembled  journal  eccentricity  and  attitude  angle  (see  Figure  2.1)  such  that  the  journal 
load  equals  the  applied  load.  Often,  a  Newton-Raphson  based  iteration  scheme  is  employed 
during  this  step  [20]. 

Once  the  bearing  is  assembled,  the  dynamic  coefficients  for  each  pad  are  determined  by  applying 
the  methods  presented  in  Section  2.3  using  the  eccentricity  and  attitude  angle  relative  to  each  pad 
in  the  pad’s  local  coordinate  system.  The  dynamic  coefficients  for  each  pad  are  computed  in  the 
local  pad  coordinate  system  (see  Figure  4.2)  assuming  the  pad  is  fixed  (i.e.,  not  free  to  rotate 
about  its  pivot  point).  However,  because  of  the  pad  rotational  degree  of  freedom,  there  is  no 
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Pad  Local  Coordinate  System  (£,n) 


Figure  4.2.  Tilting  pad  local  and  global  coordinate  systems. 

static  force  in  the  77-direction  and  the  dynamic  force  in  this  direction  is  a  function  of  the  pad’s 
rotational  effect  on  the  fluid  film  as  well  as  its  inertia  and  vibration  frequency. 


4.3  Synchronous  Reduction  of  Dynamic  Coefficients 


At  this  point,  an  assumption  of  the  pad  vibration  frequency  is  often  made  and  the  coefficients  are 
reduced  from  8+KW  to  8  (4  stiffness  and  4  damping)  coefficients  so  they  can  be  incorporated 
into  a  low-order  structural  model.  The  most  commonly  used  assumption  is  that  of  synchronous 
vibration  of  the  pads  with  the  rotational  speed  of  the  journal,  which  is  often  referred  to  as  the 
synchronous  reduction  of  the  dynamic  coefficients.  This  method  was  originally  proposed  by 
Lund  [13].  The  pertinent  equations  are  as  follows  (Additional  details  can  be  found  in  References 
11  and  13.): 


k»  =  kg  cos 2  ys  +  k'nn  sin2  ys -  (k'^  +  k'nS )cos y/ sin  iff 
kxy  =  *4  cos  V  -  k'n(  sin  V  +  {k’((  ~k'm)cosy/sinys 
kyx  =  k'„{  cos2  yf  —  k'fy  sin2 y/  +  (k'i{  -  k'nn  )cos y/ sin  iff 
kyy  =  k'nn  cos2  iff  +  k'({  sin 2  (/  +  (k^  +  k'ni)cosyr  sin  y/ 
tub ^  =  tub'cc  cos2  iff +  tub'nn  sin2  i/s-(njb'cn  +  tub'n  *)cos  ys  sin  ys 
rub ^  =  tub\v  cos 2  y/  -  tub'ni  sin 2  yr  +  (tub'g  -  tub'nn  )cos  y/ sin  y/ 
rubyx  =  tub cos 2  y/  -  tub ^  sin 2  y/  +  (tub'^  -  tub'nn  )cos  y/ sin  y/ 
wbyy  =  tub'w  cos 2  yr  +  tub ^  sin 2  y/  +  (i tub ^  +  tub'n,  )cos y/ sin  y/ 
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where 


~  kg  {pkfr  +  Qtvbfri  )knf  {ykfy  pmb^)mbn^ 
k'tn  =~Mm2(pkin  +  qmb4r>) 

Ki  =~A 1m2(pkn{+qwbni) 
k'nn  —  —M  nr2  (l  +  pM  nr2 ) 

Vb'ft  =  -  (pk{n  +  q&b'r,  )p*>n(  +  (qkfr  ~  P^bfn  Kf  ’ 

wb\n  =  Mar2(qkin  -  pmb^ ) 

wb'nn=(\im2)-q 

and 


P  = 


knr)  -Mm- 


q  = 


(k^-Mm2)2  +(mbnn)2 

_ ^nn _ 

(k„n  -Mm2)2  +(mbnn)2 


(4.4) 


The  process  of  reducing  these  coefficients  in  this  manner  introduces  a  frequency  dependence  of 
the  coefficients  as  shown  in  Equation  4.3. 


4.4  Coefficients  for  Explicitly  Modeled  Tilting  Pads 

To  avoid  the  complications  associated  with  using  reduced  coefficients  (i.e.,  frequency  dependent 
coefficients  and  no  formal  method  to  circumferentially  distribute  the  coefficients  for  a  high-order 
structural  model),  the  pads  can  be  modeled  explicitly  in  a  finite  element  model  and  the  rotor  and 
stator  can  be  coupled  through  the  pads  using  distributed  dynamic  coefficients.  The  coefficient 
distribution  is  accomplished  through  the  application  of  appropriate  integration  limits  in 
Equation  3.10.  With  this  approach,  the  pad  rotational  degrees  of  freedom  are  accounted  for  in 
the  finite  element  model  and  the  coefficients  are  not  frequency  dependent9.  Note  that  the 
coefficients  are  computed  in  a  local  coordinate  system  and  then  transformed  to  the  global 
coordinate  system  using  a  transformation  matrix,  T,  such  that  K  =  [T]T[A:][r]  and  B  =  [71T[b][r], 
where  K  and  B  are  the  coefficients  in  the  global  coordinate  system,  and  k  and  b  are  the 
coefficients  in  the  local  coordinate  system.  The  transformation  matrix,  T,  is  defined  such  that 
r  =T7  where  r  is  a  displacement  in  the  global  coordinate  system  and  r  a  displacement  in  the 
local  coordinate  system  [25].  For  the  coordinate  systems  defined  herein  (i.e.,  angles  measured 


g  Note  that  the  coefficients  are  a  function  of  the  journal  rotation  rate  because  this  has  an  effect  on  the  fluid  film 
pressure  characteristics.  However,  for  a  given  journal  speed.  Q,  all  of  the  journal  bearing  dynamic  coefficients, 
except  for  the  tilting  pad  reduced  coefficients,  are  independent  of  vibration  frequency. 
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from  the  negative  ^-direction),  the  transformation  matrix  is  computed  as  follows  for  each  pad 
coordinate  system: 


T,= 


-  cos  t//j 
sin{/( 


-  sin  y/i 

-  cos  y/i 


i  =  1  to  N. 


(4.5) 


When  computing  the  dynamic  coefficients  for  a  tilting  pad  bearing,  it  is  important  to  recognize 
that  the  final  results  are  a  function  of  the  bearing  load,  direction  of  applied  load  (e.g.,  load-on- 
pad  versus  load-between-pads),  journal  operating  speed,  as  well  as  the  other  fluid  film 
parameters  defined  above.  Additional  effects  on  the  computed  equivalent  bearing  coefficients 
such  as  those  caused  by  a  non-rigid  pivot  and  a  flexible  pad  have  not  been  addressed  here  but  are 
discussed  in  the  literature  [23,24].  Note  that  pads  for  which  there  is  no  converging  film,  called 
unloaded  pads,  generally  do  not  contribute  to  the  overall  stiffness  but  do  in  some  cases 
(depending  on  their  pivot  design)  contribute  to  the  damping  [22]. 


4.5  Application  of  Dynamic  Coefficients  to  Finite  Element  Structural  Models 

The  unsymmetric  bearing  stiffness  matrices  that  often  result  for  bearing  models  can  be  difficult 
to  implement  in  structural  finite  element  models.  The  presence  of  an  unsymmetric  matrix  causes 
the  system  to  have  complex  eigenvalues,  and  thus  a  complex  eigensolution  may  be  required. 
Standard  methods  for  incorporating  such  matrices  into  a  finite  element  code  are  available.  For 
instance,  the  direct  matrix  input  cards  (DMIG)  of  COSMIC/,  MSC/  and  CSA /  Nastran  can  be 
used  to  implement  any  type  of  mass,  stiffness,  or  damping  matrix  (symmetric  or  unsymmetric, 
real  or  complex)  using  the  M2GG,  K2GG,  or  B2GG  case  control  commands  for  symmetric 
matrices — for  all  solution  sequences — or  the  M2PP,  K2PP,  or  B2PP  case  control  commands  for 
unsymmetric  matrices — for  dynamic  response  analyses  only  (excluding  normal  modes  analyses). 
The  reader  is  referred  to  Nastran  documentation  [26]  for  descriptions  of  these  methods  and  their 
limitations.  Note  that  the  stiffness  matrix  can  be  decomposed  into  a  symmetric  matrix 
(*„=(*  +  k 1 )/  2 )  and  an  anti-symmetric  matrix  ( kanlisym  =  [k  - kT  )/2 )  such  that  k=ksym+ka„,iSym. 

The  symmetric  matrix  can  be  added  using  the  K2GG  case  control  command  while  the  anti¬ 
symmetric  matrix  can  be  added  using  the  K2PP  case  control  command.  Adding  the  matrix  in 
this  manner  and  performing  a  modal  complex  eigensolution  (Nastran  SOL  1 10)  allows  the  effect 
of  the  anti-symmetric  matrix  to  be  evaluated  in  terms  of  resonance  frequency  changes  from  the 
normal  modes  solution  to  the  complex  eigensolution.  If  the  effects  are  negligible,  then  dynamic 
results  can  be  obtained  using  the  symmetric  matrix  with  no  need  for  unsymmetric  matrix  solvers, 
which  are  generally  more  computationally  expensive. 


5  Distributed  Coefficient  Example  Problem 

An  example  problem  is  provided  to  illustrate  the  effect  of  applying  bearing  coefficients  to  a  FE 
structural  model  using:  1)  the  method  presented  herein  for  distributed  coefficients,  and  2)  a 
uniform  coefficient  distribution  to  represent  a  method  that  is  often  applied  to  such  models.  The 
example  problem  focuses  on  a  tilting  pad  bearing  design  because  of  the  complexity  associated 
with  this  type  of  bearing. 
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As  discussed  above,  the  bearing  pads  of  a  tilting  pad  bearing  are  generally  not  explicitly  modeled 
in  a  structural  model.  Instead,  the  bearing  coefficients  are  computed  using  some  sort  of 
reduction  technique  (e.g.,  synchronous  reduction)  to  account  for  the  pad  vibration  about  the 
rotational  degrees  of  freedom  and  the  coefficients  (stiffness  and  damping)  are  applied  to  the 
bearing  using  an  ad  hoc  distribution  that  may  only  capture  the  low-order  bearing  effects.  For  the 
example  problem,  the  response  of  the  model  is  computed  twice:  once  with  the  proposed  bearing 
coefficient  distribution  and  the  tilting  pads  explicitly  modeled  (the  “actual  coefficient  model”), 
and  a  second  time  using  a  uniform10  distribution  of  synchronously  reduced  coefficients  to  each 
node  on  the  rotor  circumference  (the  “reduced  coefficient  model”).  The  structure  being  modeled 
is  shown  in  Figure  5.1  with  relevant  information  shown  in  Figure  5.2.  The  outer  ring  and  tilting 
pads  are  modeled  using  steel  while  the  rotor  ring  is  modeled  using  aluminum.  The  bearing 
consists  of  eight  tilting  pads  that  are  equally  spaced  around  the  journal  circumference  and 
connected  to  the  rotor  ring  at  the  pivot  points  shown  in  Figure  5.3  with  a  pad  pivot  location  ratio, 
cdy,  of  0.5  (see  Figure  2.2).  Other  details  regarding  the  bearing  are  shown  in  Table  5.1: 


Table  5.1.  Bearing  parameters  (water  at  305  K) 


//  =  7.69x10  4  Pa  s 

Q  =  19.2  Hz 

y  =  32° 

p  -  995.  kg/m3 

L  =  0.0203  m 

o 

II 

D  =  0.323  m 

C  =2.54xl0'4m 

m  =  0.4 

II 

00 

o 

0 

Figure  5.1.  FE  model  for  the  actual  coefficient  model  showing  applied  load  and  response 

locations. 


10  The  uniform  distribution  chosen  here  has  no  physical  basis,  but  is  consistent  with  what  might  be  employed  for 
such  a  model. 
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Figure  5.2.  Example  structure  relevant  information. 


Pivot  Points 


Figure  5.3.  Tilting  pad  pivot  location. 

Because  the  tilting  pads  are  not  explicitly  modeled  when  synchronously  reduced  coefficients  are 
used,  the  model  for  the  reduced  coefficients  must  be  slightly  modified  from  what  is  shown  in 
Figure  5.1.  Two  things  should  be  noted  regarding  the  changes.  First,  because  the  pads  are 
connected  to  the  rotor  ring  in  the  model  shown  in  Figure  5.1,  concentrated  mass  elements  are 
added  to  the  reduced  coefficient  model  at  the  pivot  locations  to  account  for  the  translational 
inertia  of  the  pads  (the  rotational  inertia  about  the  pad  pivots  is  inherently  accounted  for  in  the 
synchronously  reduced  coefficients).  These  point  masses  are  shown  in  Figure  5.4.  The  rotor  and 
stator  rings  are  identical  for  both  models,  and  this  fact  enables  direct  comparisons  between  the 
results  from  each  model. 

The  first  step  in  the  analysis  is  to  compute  the  bearing  dynamic  coefficients.  Two  sets  of 
coefficients  are  computed:  one  set  of  distributed  coefficients  using  the  methods  presented  here 
with  a  distribution  to  the  FE  grids  on  the  tilting  pads,  and  a  second  set  using  the  method  of 
synchronous  reduction  and  a  uniform  distribution  to  the  FE  nodes  on  the  journal  circumference. 
In  both  cases,  the  bearing  coefficients  connect  to  the  FE  grids  on  the  inner  surface  of  the  outer 
ring.  Because  the  pad  assembly  method  is  employed  here,  both  coefficient  computations  require 
the  static  operating  characteristics  (eccentricity  and  attitude  angle)  for  each  pad  of  the  tilting  pad 
bearing. 

The  pad  static  characteristics  are  determined  by  solving  the  Reynolds  equation  and  computing 
the  hydrodynamic  force  (Equation  3.1)  for  a  variety  of  journal  eccentricities.  The  attitude  angle, 
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,  Applied  Load 


Figure  5.4.  FE  model  for  the  reduced  coefficient  model. 

0,  for  each  eccentricity  is  obtained  by  varying  0  until  the  force  in  the  77-direction  is  numerically 
zero  (the  coordinate  system  in  Figure  4.2  is  defined  such  that  the  load  is  entirely  in  the  £ 
direction).  A  sample  pressure  field  for  a  journal  eccentricity  ratio  of  0.276  and  an  attitude  angle 
of  52.2°  is  shown  in  Figure  5.5.  The  pressure  field  was  computed  using  a  finite  difference 
approximation  for  the  Reynolds  equation  with  a  discretization  of  100  in  the  ^-direction  and  60  in 
the  z-direction.  The  Reynolds  boundary  conditions  were  applied  using  an  ambient  pressure  of 
zero  for  all  boundaries.  The  net  journal  load  that  can  be  supported  by  this  pressure  field  is 
(computed  from  Equation  3.1):  =  0.371  N  and  Fn  =  -6.93xlO'10  N.  Performing  such  a 

calculation  for  several  eccentricities  yields  the  results  shown  in  Figure  5.6. 

The  pad  assembly  method,  described  in  Section  4.2,  is  then  employed  to  compute  the  overall 
journal  eccentricity  and  attitude  angle  to  support  a  desired  load  of  111.0  N.  In  doing  this,  the 
local  eccentricity  and  attitude  angle  relative  to  each  pad  is  obtained  and  the  bearing  coefficients 
can  then  be  computed  for  each  pad.  It  is  at  this  point  where  the  two  methods  diverge.  The 
synchronously  reduced  coefficients  are  computed  using  coefficients  integrated  over  the  entire 
pad  surface  and  the  methods  described  in  Section  4.3  are  employed,  while  the  actual  coefficient 
model  uses  a  distribution  of  the  coefficients  to  the  pad  FE  grids.  The  coefficients  are  distributed 
by  solving  for  the  pressure  field  over  each  pad  and  choosing  appropriate  integration  limits  in 
Equation  3.10  based  on  the  pad  grid  locations.  The  resulting  coefficients  for  all  pads  are  shown 
in  Figure  5.7  for  the  stiffness  and  Figure  5.8  for  the  damping.  As  shown  in  these  figures,  the 
stiffness  and  damping  are  largest  for  the  pad  at  180°  and  decrease  to  zero  for  the  unloaded  pad  at 
0°.  Based  on  these  results,  one  can  imagine  the  errors  that  would  be  associated  with  a  uniform 
distribution  of  these  coefficients  over  the  journal  circumference. 

Because  the  reduced  coefficients  incorporate  the  effects  of  the  pad  rotation,  an  equivalent  pad 
mass  is  required  for  the  reduced  coefficient  computation.  The  pad  mass  is  computed  as  follows 
for  the  steel  pads: 
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(5.1) 


/  _6.53xlO~bkg-m2 
K  ~  (0.165%)2 


2.50jd0‘4/:g . 


The  resulting  coefficients  after  performing  the  reduction  are  shown  in  Table  5.2.  These 
coefficients  are  equally  distributed  to  the  FE  grids  on  the  journal  circumference  of  the  reduced 
coefficient  model. 
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Figure  5.6.  Sample  fixed  pad  static  characteristics  (W  =  journal  load). 
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Table  5.2.  Synchronously  reduced  coefficients. 


kxx  =  9.71xlOb  N/m 

kxv  =  -77.8  N/m 

kvx=  132.  N/m 

kvv  =  3.36xl05  N/m 

bxx=  1.98x1 04N/(nvs) 

bxv  =  0.797  N/(m  s) 

bvx  =  -1.19  N/(m-s) 

bvv  =  1. 54x10'*  N/(m-s) 

It  is  interesting  to  note  that  it  is  the  inclusion  of  the  pad  equivalent  mass  in  this  calculation  that 
causes  the  cross  coefficients  to  become  non-zero.  Performing  the  same  computation  with  no  pad 
mass  yields  zero  cross  coefficients  as  shown  in  Table  5.3.  The  pad  mass  does  have  an  effect  on 
all  of  the  coefficients  in  these  tables,  but  the  differences  are  small  relative  to  the  diagonal  terms 

(i.e.,  kxx,  kyy,  bxx,  and  byy). 
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Figure  5.8.  Distributed  damping  coefficients  for  the  actual  coefficient  model. 


Table  5.3.  Synchronously  reduced  coefficients  using  no  pad  mass. 


kxx  =  9.71xlObN/m 

o 

II 

>■ 

X 

O 

II 

X 

>• 

kvv  =  3.36xl05  N/m 

bxx  =  1.98xl04  N/(m-s) 

bXy  =  0 

bVx  =  0 

byy  =  l^xlO'1  N/(ms) 

Also  note  that  the  pad  located  at  i//=  0  is  unloaded  for  this  bearing  and  thus  does  not  contribute 
to  the  bearing  coefficients  (but  the  damping  could  contribute  depending  on  the  pivot  design).  If 
this  pad  were  fully  loaded  (i.e.,  a  converging  film  over  the  entire  surface),  its  contribution  would 
cause  the  direct  coefficients  in  Table  5.2  and  Table  5.3  to  be  equal  because  of  the  bearing 
symmetry  (i.e.,  k^=kyy,  and  bxx=byy). 

Accelerance  (acceleration/force)  transfer  function  plots  for  a  radial  excitation  on  the  outer  ring 
and  a  response  location  on  the  rotor  (see  Figure  5.1  and  Figure  5.4)  are  shown  in  Figure  5.9  for 
both  models.  These  results  were  obtained  by  performing  a  frequency  response  analysis  using 
MSC/Nastran  with  the  bearing  coefficients  added  via  DMIG  cards  [26].  As  can  be  seen  in  this 
figure,  the  transfer  accelerances  have  similar  characteristics  but  exhibit  differences  in  both 
amplitude  and  resonance  content  for  the  two  models.  These  differences  reflect  the  method  by 
which  the  bearing  coefficients  are  applied  since  the  models  are  otherwise  identical.  The  mean 
accelerance  (response  points  averaged  over  the  rotor  inner  surface)  and  corresponding 
circumferential  Fourier  decomposition  for  each  model  due  to  the  same  drive  point  as  before  is 
shown  in  Figure  5.10.  These  plots  show  a  significant  amount  of  cross  modal  coupling  for  the 
actual  coefficient  model  relative  to  the  reduced  coefficient  model  for  the  resonances  near  300 
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Figure  5.9.  Comparison  of  transfer  accelerance  magnitude  using  the  proposed  coefficient 
distribution  and  a  uniform  distribution  of  synchronously  reduced  coefficients. 

and  800  Hz.  The  Fourier  decomposed  results  are  obtained  using  the  discrete  Fourier  transform 
(DFT)  applied  along  the  circumferential  direction  for  each  frequency  as  follows  [27]: 


1  NNS  1  N- 1  -i(— 


w 


NNS  N 


n=0 


(5.2) 


where  a\v\  is  the  average  accelerance  at  a  specific  frequency  for  the  Fourier  component  v,  NNS 
is  the  number  of  nodes  per  segment  of  the  structure  (the  structure  is  divided  into  N  equal 
segments  along  the  circumferential  direction),  and  a\i,fi\  is  the  accelerance  at  node  i  of  segment 

ju.  Because  only  the  first  N/2  Fourier  components  can  be  resolved  for  an  /V-point  DFT  and  these 
coefficients  are  “mirrored”  about  N/2,  the  mirrored  coefficients  are  summed  together  such  that 

a[l]  =  abs(a[l])  +  abs(a[/V  - 1]),  a\l\  =  abs(a[2])  +  abs(a[/V  -  2]),  etc.  to  ensure  the 
decomposed  results  sum  together  to  yield  the  overall  response.  Note  that  if  N  is  an  even  number, 
the  a[N/2 ]  component  is  not  mirrored.  Here,  N  is  chosen  to  be  eight  based  on  the  structure 
having  eight  times  symmetry,  which  is  a  result  of  having  eight  tilting  pads  equally  spaced  around 
the  rotor  circumference.  This  choice  of  N  for  the  reduced  coefficient  model  is  sufficient  to 
resolve  all  of  the  Fourier  components  because  the  structure  is  axi-symmetric  with  eight  times 
symmetry.  The  actual  coefficient  model  is  not  axi-symmetric  due  to  the  presence  of  a  non- 
uniform  bearing  coefficient  distribution.  Thus,  a  larger  value  of  N  should  be  chosen  for  the 
actual  coefficient  model  to  decompose  the  response  with  minimal  aliasing  effects.  However,  the 
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choice  of  N=  8  is  sufficient  for  both  models  to  show  the  effect  of  the  bearing  distribution  method 
on  the  Fourier  content1 1  (see  Figure  5.10). 


Actual  Coefficient  Model 


8 

c 


O  Reduced  Coefficient  Model 


Frequency  (Hz) 


-  total 

N=0 

N=1 

-  N=2 

N=3 
.  N=4 


Figure  5.10.  Mean  rotor  accelerance  decomposition  into  circumferential  Fourier 

components. 


A  summary  of  resonance  frequencies  and  structural  damping  loss  factors  for  the  two  models  is 
provided  in  Table  5.4  for  resonances  up  to  1  kHz.  Several  aspects  of  these  results  should  be 
noted.  First,  numerous  modes  are  overdamped  due  to  the  presence  of  the  bearing  damping,  and 
thus  do  not  appear  in  this  table.  Second,  several  modes  have  a  damping  consistent  with  the 
assumed  material  damping  of  0.001  (e.g.,  298.,  429.,  and  818.  Hz  for  the  actual  coefficient 
model,  and  310.,  429.,  and  874.  Hz  for  the  reduced  coefficient  model).  These  are  torsional 
modes,  indicated  by  a  "T"  in  Table  5.4,  that  are  not  affected  by  the  bearing  coefficients  because 
the  coefficients  act  only  in  the  radial  direction  and  the  centerline  rotor  and  outer  ring  grids  where 
the  bearing  coefficients  are  applied  have  displacements  only  in  the  axial  direction  for  these 


1 1  An  evaluation  using  various  values  for  N  for  the  actual  coefficient  model  decomposition  has  shown  that  aliasing 
effects  have  a  negligible  impact  on  the  results  reported  in  Figure  5.10. 
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modes.  The  differences  between  the  resonance  frequencies  of  each  model  for  these  modes  are  a 
result  of  neglecting  the  pad  rotational  inertia  about  the  r  and  6  directions  in  the  reduced 
coefficient  model.  Third,  modes  for  which  the  rotor  is  spatially  comprised  of  a  single 
circumferential  Fourier  component  are  identified  in  this  table  with  the  variable  "N"  (bending 
modes  only).  Modes  without  such  a  designation  are  largely  affected  by  the  non-uniform 
coefficient  distribution  and  are  thus  comprised  of  multiple  Fourier  components,  which  are  not 
shown  in  the  table.  Next,  resonances  corresponding  to  the  peaks  in  Figure  5.9  are  shown  in  bold 
face  font  in  Table  5.4.  For  these  modes,  both  the  rotor  and  outer  ring  are  comprised  of  the  same 
Fourier  component  shown  in  the  table.  Moreover,  it  is  apparent  from  the  damping  in  this  table 
that  the  peak  height  near  800  Hz  in  Figure  5.9  is  not  due  to  disparate  damping  levels,  but  due  to 
structural  dissimilarities  (the  transfer  accelerance  peak  and  loss  factor  for  the  reduced  coefficient 
model  are  both  higher  than  those  for  the  actual  coefficient  model).  Finally,  numerous  repeated 

Table  5.4.  Resonance  frequencies  and  damping  levels. 


Actual  Coefficient  Model 

Reduced  Coefficient  Model 

Resonance 
Frequency  (Hz) 

Loss  Factor 

Resonance 
Frequency  (Hz) 

Loss  Factor 

3.72  (N=l) 

16. 

7.72  (N=l) 

6.0 

19.5  (N=l) 

6.7 

153. 

0.62 

169.  (N=2,x2) 

2.1 

261.  (N=2) 

0.068 

267.  (N=2) 

0.73 

288.  (N=2) 

0.022 

286.  (N=2,x2) 

0.018 

287.  (N=2,x2) 

0.016 

298.  (T,x2) 

0.001 

310.  (T,x2) 

0.001 

298.  (N=2) 

0.026 

299.  (N=2) 

0.022 

429.  (T,x2) 
(outer  ring) 

0.001 

429.  (T,x2) 

(outer  ring) 

0.001 

548. 

0.48 

615. 

0.54 

696.  (N=3) 

0.18 

746.  (N=3) 

0.27 

747. 

0.62 

788.  (N=3) 

0.031 

809.  (N=3) 

0.025 

794.  (N=3) 

0.044 

798.  (N=3) 

0.041 

816. 

0.053 

818.  (T,x2) 

0.001 

874.  (T,x2) 

0.001 

913. 

0.42 

976. 

0.22 
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resonances  occur  for  the  reduced  coefficient  model  due  to  the  model’s  uniform  bearing 
coefficient  distribution  and  the  structure’s  axi-symmetry  (these  are  indicated  by  a  "x2"  next  to 
the  resonance  frequencies  in  Table  5.4).  Repeated  resonances  are  also  present  in  the  actual 
coefficient  model,  but  only  for  modes  not  affected  by  the  non-uniformly  distributed  bearing 
coefficients  (i.e.,  the  torsional  modes  described  previously). 


6  Summary  and  Conclusions 

A  review  of  methods  currently  employed  for  modeling  static  and  dynamic  characteristics  of  fluid 
film  bearings  has  been  provided.  Pertinent  literature  has  been  reviewed  and  summarized  with  an 
effort  placed  on  highlighting  the  nuances  associated  with  such  models.  In  nearly  all  cases,  the 
literature  discusses  the  use  of  dynamic  coefficients  for  rotor-dynamics  models,  in  which  a  single 
stiffness  and  damping  matrix  is  used  to  connect  an  individual  rotor  node  to  an  individual 
stator/bearing  node.  The  focus  of  this  work  was  on  developing  dynamic  coefficients  for  a 
structural  finite  element  model,  which  uses  a  distribution  of  the  coefficients  over  the  journal 
circumference.  The  proposed  method  is  applicable  to  all  fluid  film  bearings.  The  method  is 
approximate  in  the  sense  that  the  effect  of  a  flexible  journal  and/or  bearing  is  not  incorporated  in 
the  resulting  coefficients. 

For  the  case  of  a  tilting  pad  bearing,  rotor-dynamics  models  often  use  a  reduced  set  of 
coefficients  such  that  the  individual  pads  need  not  be  modeled.  The  process  of  reducing  the 
coefficients  causes  them  to  become  frequency  dependent  since  a  frequency  for  the  pad  rotational 
degrees  of  freedom  must  be  assumed.  As  discussed  herein,  by  explicitly  modeling  the  pads  in, 
for  example,  a  structural  finite  element  model,  the  pad  rotational  degrees  of  freedom  are 
accounted  for  in  the  finite  element  model  and  thus  the  bearing  coefficients  are  no  longer 
frequency  dependent. 

A  brief  discussion  regarding  the  application  of  the  computed  dynamic  coefficients  for  a 
structural  finite  element  model  is  also  provided.  Because  of  the  unsymmetric  nature  of  the 
stiffness  coefficients  (the  damping  coefficient  matrix  should  always  be  symmetric),  there  may  be 
a  need  to  directly  input  the  matrices  to  the  finite  element  model,  depending  on  the  available 
elements  for  the  solver  being  used  (e.g.,  ABAQUS,  ANSYS,  or  NASTRAN).  It  was  also 
suggested  that  symmetric  and  anti-symmetric  matrices  be  created  from  the  unsymmetric  matrices 
to  facilitate  the  solution  process. 

The  importance  of  appropriately  distributing  the  computed  bearing  coefficients  is  demonstrated 
through  a  numerical  study  wherein  a  comparison  of  results  for  a  tilting  pad  bearing  model  are 
computed  for  two  different  implementations  of  the  bearing  coefficients.  The  two  coefficient 
implementations  are:  1)  a  synchronous  reduction  of  the  coefficients  and  a  uniform  distribution  to 
the  journal  circumference  in  a  model  with  no  tilting  pads,  and  2)  a  model  using  the  distribution 
method  proposed  herein  which  incorporates  explicit  structural  models  of  the  tilting  pads  and  a 
distribution  of  the  coefficients  to  each  pad  grid  in  the  FE  model.  A  comparison  of  results  for 
these  models  shows  substantial  differences  in  the  transfer  accelerance  (both  in  amplitude  and 
circumferential  makeup),  resonance  frequencies,  and  damping  levels. 
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Appendix:  Simplified  Solutions  to  Reynolds  Equation 


Simplified  solutions  to  Equation  2.1  are  available  for  steady-state  conditions 


=  0)  for  the 


case  of  an  infinitely  long  plain  journal  bearing  (i.e.,  ^  » 


1  =►  ^  )  and  for  the  case  of  a  short 


plain  journal  bearing  (i.e.,  c  1  =*  k<  ^  )•  The  solution  to  the  long  bearing  approximation, 
called  the  long  bearing  solution,  or  the  Sommerfeld  solution,  is  shown  in  Equation  A.l  [4]. 


P  = 


pRQ 


6e  sin/(2  +  £'cos?') 


c~  [  (2 +£2\[u72 


(A.l) 


where  C  is  the  radial  clearance,  e  is  the  ratio  of  the  eccentricity  and  the  radial  clearance  {e!C), 
0  is  the  angle  between  the  load  direction  and  the  film  rupture  location  (see  Figure  2.1),  and  pa 
accounts  for  any  supply  pressure  at  y=  G-<p=  0.  Note  that  film  rupture  occurs  at  y=  n  and  the 
active  film  falls  within  0<  y<n.  Equation  A.l  predicts  negative  pressures  for  n<  y<  0,  and 
thus  the  pressure  is  assumed  to  be  pa  within  this  range.  A  relationship  between  the  load  and 
eccentricity  is  provided  in  Equation  A. 2. 


_  fjR^QL  12 ne 
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where  W  is  the  load. 


The  short  bearing  solution,  or  the  Ocvirk  solution,  is  shown  in  Equation  A.3,  with  the  load- 
eccentricity  relationship  in  Equation  A.4  [4]. 
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where  L  is  the  bearing  width,  and  -  —  <  z  <  — . 
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The  film  rupture  location  is  determined  using  Equation  A.5. 


<p  =  tan 
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